#include "stdafx.h"
#include "Random.h"

namespace Cvekas {

Random::Random()
{
   init((uint)time(0));
}

Random::Random(uint seed)
{
   init(seed);
}

uint Random::randUInt()
{
   if(index >= n)
      update_state();

   uint result = state[index++];
   
   // Temper
   result ^= (result >> 11);
   result ^= (result << 7) & 0x9d2c5680UL;
   result ^= (result << 15) & 0xefc60000UL;
   result ^= (result >> 18);

   return result;
}

uint Random::randUInt(uint low, uint high)
{
   assert(low < high);

   return low + randUInt() % (high - low);
}

int Random::randInt()
{
   return int(randUInt());
}

int Random::randInt(int low, int high)
{
   assert(low < high);

   return low + randUInt() % (high - low);
}

float Random::randFloat(float low, float high)
{
   assert(low < high);

   // Random float in range [0, 1)
   float r = float(randUInt()) / (float(MAX_UINT32) + 1.0f);
   
   return low + r * (high - low);
}

// Box-muller transform based normal distribution
float Random::normal()
{
   static float cached = infinity;
   float result, r, t;

   if(cached != infinity)
   {
      result = cached;
      cached = infinity;
      return result;
   }

   r = sqrt(-2.0f * log(randFloat()));
   t = 2.0f * pi * randFloat();
   result = r * cos(t);
   cached = r * sin(t);
   return result;
}

float Random::normal(float mean, float variance)
{
   assert(variance > 0.0f);

   return sqrt(variance) * normal() + mean;
}

void Random::init(uint seed)
{
   state[0] = seed;
   for(index = 1; index < n; ++index)
   {
      uint t = (state[index-1] ^ (state[index-1] >> 30));
      state[index] = 1812433253UL * (t + index);
   }
}

// TODO: investigate if it's possible to do the same with 1 loop
void Random::update_state()
{
   int kk, y;
   
   const uint upper = 1UL << 31;
   const uint lower = ~upper;
   const uint matrix = 0x9908b0dfUL;

   for(kk = 0; kk < n - m; ++kk)
   {
      y = (state[kk] & upper) | (state[kk+1] & lower);
      state[kk] = state[kk + m] ^ (y >> 1); 
      state[kk] ^= (y & 1) ? 0 : matrix;
   }

   for(; kk < n - 1; ++kk)
   {
      y = (state[kk] & upper) | (state[kk+1] & lower);
      state[kk] = state[kk + (m - n)] ^ (y >> 1);
      state[kk] ^= (y & 1) ? 0 : matrix;
   }
   
   y = (state[kk] & upper) | (state[0] & lower);
   state[kk] = state[kk + (m - n)] ^ (y >> 1);
   state[kk] ^= (y & 1) ? 0 : matrix;

   index = 0;
}

} // namespace